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Abstract 

In analogy with the quantum field theory of free bosons a sim¬ 
ple integral representation is derived for recently proposed corrections 
describing the Bose Einstein effect. The saddle point approximation 
to these integrals results in a compact expression which sums effec¬ 
tively the original n! terms with an accuracy better than 2% for n > 7 
strongly correlated bosons. 
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1 Introduction 


The very nature of the Bose-Einstein correlations makes them rather difficult 
to include in the Monte Carlo event generators and consequently the problem 
has been the area of intensive studies for many years |T| -||. Recently Bialas 
and Krzywicki have proposed a simple recipe to effectively include Bose- 
Einstein correlations in the event generating Monte Carlo programs j5|. The 
basic step in their procedure consists of calculating a positive weight for each 
generated event 


n 

W n (p 1, ...,p n ) = En A iQ n (i) =perm(A), (1) 

Qn i 

where A is a correlation matrix which depends on particle momenta and on 
the choice of variables parameterizing the whole phenomenon. One com¬ 
monly used form reads 


A ik = exp 


^ ( Pi-Pk ) 2 j 


( 2 ) 


The sum in Eq.(|T]) extends over all permutations, Q n , of n elements (1,2,... ,n). 
Hence the numerical cost to calculate above weights grows like n\ and this 
limits practical applications. 

In this letter we derive a simple 2n -dimensional integral representation for 
a general permanent of a correlation matrix, perm(A). Our representation 
has two advantages. First, the integrals can be done analytically, in the sad¬ 
dle point approximation, providing relatively simple and accurate expression 
for the sum (]!]). Secondly, the integral representation has a straightforward 
probabilistic interpretation. This indicates that implementation of the cor¬ 
rection scheme proposed in Ref.[5| might be possible stochastically, and may 
advance the novel generation of MC programs. 

The motivation for this work came from the simple observation that even 
the free quantum field theory of bosons automatically incorporates sym- 
metrization. It is therefore plausible that the simple integral representation, 
analogous to that for the propagation of n identical bosons, should exist for 
the sum (|l|). In the next Setion we prove that this is indeed the case. 
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2 Integral representation 

Let us recall the Wick theorem for finite, say n, number of gaussian variables 
(0i, 02,..., 0 n ) = 0 T . If 0 is distributed according to [] 

P(0) ~ exp ^0 T G^), (3) 

then moments of 0’s are given by the sum over all contractions 

^ 0*10*2 ' ' ‘ 0*fe (0*10*2 0*30*4 ' ' •) d" (4) 

.-- ' 

(0*1 0*20*3 0*4 - - •) 4“ • - - ; 
where a single contractions reads 


Mt= (<?-%, 


( 5 ) 


for a non-singular matrix G'. 

To match the combinatorics of contractions with that of simple permu¬ 
tations required in Eq.(!l[) we introduce 2 n gaussian variables 
(0i,..., 0 n ), (0i,..., 0 n ) = (0 T , 0 T ) = T 7 distributed accordingly to 


p(vp) ~ exp 




with the 2n x 2n matrix 


M = 


0 

A -1 



( 6 ) 

( 7 ) 


and A as given by Eq.(|2|). Since, by construction, 0 and 0 variables are not 
coupled within themselves, all non-zero contractions occur between 0’s and 
0’s only. Hence it follows from the Wick theorem (|5|) that 

n 

< 11(0^*) >= W n (pi,...,p n ). (8) 

i 

1 Where evident the vector label "* will be omitted. 
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since M 1 = 



This is the integral representation sought of in 


the Introduction. However the choice of variables, although simple concep¬ 
tually, is not very practical, since the form is not positive definite, 

and consequently the paths of fl" dijj t d(f) l integrations should lie in the com¬ 
plex plane of the corresponding variables. This deficiency is removed and 
the integrand further simplified by the following series of transformations of 
variables (j, k — 1,... n) 


‘A j 


Ah 


1 / \ 
-/=( Uj + IVj), 

'U'j OjkCk i 

Ci = 

yfcxj, 

( 9 ) 

-j={uj — ivj ), 

Vj OjkT)ki 

Vj = 


( 10 ) 


where the orthogonal matrix O transforms A into the diagonal form 0 T AO = 
Diag[Xi ,..., A n ]. One obtains 


w n = /“exp + n /(^e ( '‘ ) ) 2 + fa'e < " , ) 2 )n^. (11) 

Vectors are given by 

e ^ = O^yfXi, (12) 

and all integrals are along the real axis. It turns out that all eigenvalues 
of the correlation matrix, Eq.(^), are nonnegative. The more general case 
of arbitrary eigenvalues can be readily dealt with by the appropriate Wick 
rotations which will be different for different coordinates. 

Equation (|Tl ) represents the sum of n\ terms in ([!]) by the 2n-fold inte¬ 
gral of the gaussian type. The whole dependence on the original momenta 
(pi,... ,p n ) is contained in vectors which, in view of Eq.(|l^), are closely 
related to the eigenvectors of the correlation matrix A. 

The integral (TO) can be calculated by the Monte Carlo technique. We 
have found, however, that the importance sampling according, to the normal 
distribution, is not sufficient to produce a decent estimate. The product 
of n factors dramatically modifies the relatively mild gaussian dependence 
and should be included in the weight. We have therefore solved the integral 
analytically in the saddle point approximation. This gives rather satisfactory 
results as discussed in the next Section. 
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3 Saddle point approximation 


3.1 The action and the saddle point equations 


The interplay between the gaussial fall off and the polynomial rise of the two 
factors in the integrand of Eq. (|lTD , produces a maximum far from the origin. 
This maximum becomes more and more narrow with increasing n. Thus, we 
expect that the saddle point approximation will be better for higher n, just 
where the exact summation 0 fails in practical terms. 

If one writes the integrand in Eq. (|TT|) as 

1 

I (x,y) = j^-exp(-S(x,y)), (13) 

then the ’’action” S(x,y ) reads 

1 n 

S{x,y) = -(x 2 + f) - J2ln (xl + yl), (14) 

where the greek subscript denotes the scalar product z fJi = (z ■ e ^), z = 
x , y. Then the location of the saddle point (or points) is determined by the 
following system of nonlinear equations 


x 

y 


2E 

n 

2 E 


pit 1 ). 


(15) 


+ vl) ’ 

pit 1 ). 

Vn 

(16) 


+ vl)' 


In general many solutions may exist and there is no simple method to find 
all of them. The problem is partially alleviated by identifying the continuous 
symmetries of the integrand. This is also necessary for the proper application 
of the saddle point technique. Therefore we will first discuss the symmetries 
of the problem and identify corresponding zero mode. 


3.2 Global U(l) symmetry and gauge fixing 

The integrand in Eq. ([TT|) is invariant under the ’’horizontal” U(l) rotation 

x(o>) = x cos (a) + y sin (a), (17) 

y(ct) = — ifsin (a) + ycos (a), (18) 
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with the same a for every coordinate x t , y*. Consequently there is no unique 


solution to Eqs.(15-16) and one has to fix this freedom before the saddle 
point technique can be used. We follow the standard procedure known from 
the field theory. In this language the symmetry (|T8|) is the global one, and 
the single ’’gauge fixing” condition is required. We choose the symmetric 
gauge 

£ 2 = f, (19) 

and insert the following representation of unity under the integral (|TT|) 


r2n 

A (x,y) / 5(x(a) 2 — y(a) 2 )da = 1. 

Jo 


The gauge invariant function A (x,y) follows from Eq. 

A(x,y) = l(f 2 + ?). 


( 20 ) 


( 21 ) 


3.3 The formula 


With the aid of Eqs. © and © the volume of the gauge orbit (27r) can be 
readily factored out and the zero mode corresponding to the symmetry (|T8|) 
becomes fixed. Further calculations are straightforward and the final result 
reads 

W n = 2 1_n |s 0 |jjy- exp(-S(s 0 ,s 0 )), (22) 

V n, K 

where the product is taken over all, but one (equal to zero), eigenvalues of 
the 2n dimensional covariance matrix. 


( 2 ) _ d 2 S(x,y) l 
jk dzjdzk l(so ’ So) 


(23) 


evaluated at the saddle point (s 0 ,So)- For completness we quote the explicit 
form of S^ for arbitrary x and y. 


S (2) (x,y) 


S( xx ) \ 

S(v x ) s (w) / ’ 


(24) 
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q(xx) 
^ik 


c( x v) _ c(v x ) 
^ik ^ik 


c(yy) 


71 n /2 _ /- y »2 

S lt -2j2efef y “ * 


(■ X l + y 2 a ) 2 ’ 


x v' (ao (ao 
4 2^ e i e * 


**-2£e} 


+ l/ 2 ) 2 ’ 

~2 _ 2 
(aO _(m) x a» 


(^m + 2/m) 


2~\2 ' 


(25) 

(26) 
(27) 


Condition © does not entirely fix the freedom (|T8|). For each symmetric 
solution, (s, s) say, there exist also three rotated solutions (— s, s) (s, —s) 
and (— s, —s) which belong to the orbit ( ]l8l) and satisfy Of course all 

four (equal) contributions are included in the final result. Chosing different 
gauge may lead to the different discrete degeneracies, but the final result 
is not changed. More important, we have found that, except for the above 
degeneracy, there is no other, essentially different, solution of the saddle 
point equations. This seems to be a general property of our S(x,y ) which 
was checked for a number of choices of correlated momenta (p\ ,... ,p n ). 

Finally, a few comments about the numerical solution of the system (|H>| - 
B- Since our experiments indicate that after a complete gauge fixing there is 
only a single maximum, any numerical routine should perform satisfactorily. 
We have used a simple iterations of Eqs.(|F|jT6]) which quickly converged to 
the cycle between two parallel vectors. Then, a repeated mixture of the 
arithmetic average between the two followed by one iteration of the map (|T5|- 


converged rapidly to the unique solution. Condition (pT|) was assured by 
the choice of the symmetric starting point. 


4 Results and further applications 

Table I contains comparison between the saddle point formula, Eq.(^), and 
the exact one (]]]). As expected, the accuracy of the saddle point calculation 
increases with increasing n reaching the level of 2% at n = 8 — 9. Con¬ 
sequently, both methods applied complementarily would give a satisfactory 
approach for all n. 

Results presented in Table I were obtained for the one dimensional mo¬ 
menta generated randomly from the interval 0 < pk < 2 and with the width 
in Eq. @ a — 1/\/2. This corrresponds to very strongly correlated particles. 
In real applications, few particles are clustered within a range of typically 
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measured er ~ O.lGeVh This situation corresponds to the block diagonal 
(after suitable permutation) form of matrix A - each block describing one 
cluster of identical bosons, ft can be readily seen that, for exactly block 
diagonal (reducible) A, the integral ([□]) factorizes 

W n =W n W n ,...W nNM , r . (28) 

This factorization follows even more clearly from the definition (Jf). On the 
other hand the saddle point formula applied directly to the whole matrix A 
would fail since it does not accomodate multiple zero modes (one for each 
cluster) which occur in the reducible case. The remedy of this difficulty is 
evident. One should first perform reduction of a given event into irreducible 
clusters and then, depending on the size of the cluster, apply either exact or 
saddle point formula to each block separately. The reduction procedure will 
contain an essential parameter, say e, which can be defined as a tolerance 
beyond which bosons would be considered uncorrclated. This tolerance can¬ 
not be too small since then a quasi block diagonal matrix would be qualified 
as a single irreducible cluster. Therefore the value of epsilon is determined 
by the performance of the saddle point formula for weakly broken additional 
U(l) symmetries. Our numerical experiments indicate that the present ’’one 
loop” formula works for e > 0.05 — 0.1. 

Further improvement in the accuracy of the saddle point result can be 
easily achieved by calculating higher corrections. Including for example the 
two additional terms of the Taylor expansion around the saddle should not 
increase much the computational effort but would reduce the errors consider¬ 
ably. More important, this should allow to decrease the tolerance parameter 
e discussed above. 

Another path for an improvement is the Monte Carlo calculation of the 
corrections to the saddle point formula. As mentioned Section 2, MC calcula¬ 
tion of the average HD was not satisfactory, because the samples were drawn 
from the wrong distribution (centered gaussians). However one can use the 
saddle point approach to locate the maximum of the integrand and then use 
the saddle point approximation of the integrand as the proper distribution 
of generated samples. In that case the observable (the ratio between exact 
and saddle point integrands) will be a slowly varying function of x and y and 
statistical estimates will be reliable. 

Finally, there exists an intriguing possibility of defining a new family of 
event generating MC programs which would have the corrections (|lj) included 


on the stochastic basis. That may be possible since the integrand in our rep¬ 
resentation © is positive. Therefore it might be feasible to modify existing 
event generators by incorporating that factor. 

To summarize, we have derived an exact integral representation for the 
Bose-Einstein weights, which until now have been computed by summing 
over n\ terms. Resulting 277,-dimensional integral was calculated by the sad¬ 
dle point technique, and obtained formula approximates the exact sum to 
less than 2% for more than 7 correlated bosons. Various applications and 
improvements have been outlined. 

I would like to thank K. Fialkowski for the discussion. This work is 
supported in part by the Polish Committee for Scientific Research under the 
grant no. 2P03B19609. 
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n 


t saddle (^6c) 

TJ /'exact 
vv n 

t exact (^6c) 

error (%) 

a 


1.000 

0.4 

1.125 


12.5 

0.71 


1.878 

0.6 

2.094 


10.3 

0.68 

5 

13.144 

1.2 

14.196 

0.1 

7.4 

0.64 

7 

649.01 

2.0 

666.78 

5.0 

2.7 

0.74 

8 

8628.2 

2.5 

8748.2 

44.7 

1.3 

0.82 

9 

43688. 

3.2 

44391. 

443. 

1.5 

0.79 




n! 




10 

5.223e04 

3.9 

3.63e06 



0.65 

15 

2.216e09 

7.0 

1.31el2 



0.65 

20 

1.034el5 

11.4 

2.42el8 



0.68 

25 

2.557e22 

16.4 

1.55e25 



0.77 

30 

3.118e26 

22.5 

2.65e32 



0.63 

35 

4.112e32 

30.1 

1.03e40 



0.61 

40 

3.449e39 

36.0 

8.16e47 



0.62 


Table 1: Saddle point results and the exact summation for different number 
n of identical, strongly correlated, bosons. Timing (in seconds) of the Math- 
ematica runs on a HP735 workstation are also quoted. The last column gives 
the result of a simple fit W n = n\a n , where the parameter a 2 is a measure of 
the average value of the matrix elements of A. 
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